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Shot  peening  induced  compressive  residual  stresses  are  often  introduced  in  Ni  base  superalloy 
components  to  help  prevent  or  retard  surface  fatigue  crack  initiation  and  early  growth  at  near  surface 
inclusions.  In  certain  cases  these  compressive  residual  stresses  can  shift  the  fatigue  crack  initiation  site 
from  surface  to  sub  surface  locations.  However,  the  ability  to  computationally  predict  the  improvement 
in  fatigue  life  response  and  scatter  due  to  induced  compressive  residual  stresses  are  lightly  treated  in  the 
literature.  To  address  this  issue,  a  method  to  incorporate  shot  peened  residual  stresses  within  a  3D 
polycrystalline  microstructure  is  introduced  in  this  work.  These  residual  stresses  are  induced  by  a 
distribution  of  fictitious  or  quasi  thermal  expansion  eigenstrain  as  a  function  of  depth  from  the 
specimen  surface.  Two  different  material  models  are  used,  a  J2  plasticity  and  a  crystal  plasticity  model. 
First,  the  J2  plasticity  model  with  combined  isotropic  and  kinematic  hardening  is  used  to  determine  the 
distribution  of  quasi  thermal  expansion  eigenstrain  as  a  function  of  depth  from  the  surface  necessary  to 
induce  the  target  residual  stress  profile  within  the  microstructure.  This  distribution  of  quasi  thermal 
expansion  eigenstrain  is  then  used  within  a  crystal  plasticity  framework  to  model  the  effect  of 
microstructure  heterogeneity  on  the  variability  in  residual  stresses  among  multiple  instantiations.  This 
model  is  verified  with  experimental  X  ray  diffraction  (XRD)  data  for  scatter  in  residual  stresses  for  both 
the  initial  microstructure  and  after  a  single  load/unload  cycle. 

Published  by  Elsevier  Ltd. 


1.  Introduction 

3.3.  Motivation  for  modeling  residual  stresses 

The  beneficial  effects  of  compressive  surface  residual  stresses 
on  high  cycle  fatigue  (HCF)  response  have  been  well  documented 
in  the  literature  [1  4],  For  Ni  base  superalloys,  a  shift  from 
surface  dominated  to  subsurface  dominated  fatigue  crack  initia 
tion  sites  exists  for  the  transition  from  low  cycle  fatigue  (LCF)  to 
HCF  regimes  [5  9],  A  similar  surface  to  subsurface  transition  has 
been  reported  for  titanium  alloys  [10,11]  and  high  strength  steels 
[12,13],  In  the  transition  fatigue  regime  (cycles  to  failure,  Nf~  lx 
104  5  x  105  cycles),  subsurface  initiated  fatigue  cracks  tend  to 
require  more  cycles  to  failure  as  compared  to  surface  initiated 
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fatigue  cracks  [5],  Accordingly,  surface  compressive  residual  stres 
ses  are  often  introduced  in  Ni  base  superalloy  components  to  help 
prolong/retard  fatigue  crack  initiation  [14]  and  early  growth  at 
near  surface  inclusions  and  potentially  shift  fatigue  crack  initiation 
sites  from  surface  to  subsurface  locations  [1,2  to  increase  fatigue 
life.  However,  compressive  residual  stresses  are  usually  only  useful 
in  the  transition  fatigue  and  HCF  (Nf  >  5  x  10s  cycles)  regimes 
since  residual  stress  relaxation  at  higher  applied  stress/strain 
values  can  eliminate  the  effectiveness  of  compressive  residual 
stresses  on  fatigue  life  [14  17], 

Compressive  surface  residual  stresses  can  be  applied  via  multi 
pie  techniques  (shot/gravity  peening,  low  plasticity  burnishing, 
laser  shock  peening,  etc.  [18]).  Shot  peening  is  the  most  commonly 
used  technique  in  industry  to  induce  compressive  residual  stresses 
at  the  surface,  and  is  the  focus  of  this  residual  stress  study.  During 
the  shot  peening  process  multiple  high  velocity  shot  beads  impact 
the  surface  forming  multiple  indentations  and  inhomogeneous 
compressive  elastic/plastic  deformation  of  the  near  surface  layer. 
As  such,  the  resulting  residual  stress  profile  is  dictated  by  (1)  the 
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interaction  between  multiple  indentations,  (2)  localized  con 
strained  compression,  (3)  strain  rate  sensitivity  and  elastic/plastic 
response  of  the  material,  and  (4)  local  microstructure.  In  modeling, 
these  discrete  impingement  events  can  either  be  modeled  expli 
citly  or  idealized  as  a  collective  elastic/plastic  constrained  com 
pressive  loading/unloading  event  as  schematically  shown  in  Fig.  1. 
The  idealization  of  the  multiple  impingement  shot  peening  pro 
cess  as  a  single  deformation  event  is  adopted  in  this  paper. 

While  the  influence  of  residual  stress  on  fatigue  life  has  been 
reported  extensively  in  the  literature,  the  ability  to  computation 
ally  predict  the  improvement  in  fatigue  resistance  and  change  of 
scatter  due  to  induced  compressive  residual  stresses  are  lacking  in 
the  literature.  A  significant  effort  in  modeling  inclusion  and 
residual  stress  effects  in  shot  peened  martensitic  gear  steels  was 
undertaken  by  Prasannavenkatesan  et  al.  [19  22],  They  considered 
separate  3D  finite  element  method  (FEM)  models  at  discrete 
depths  subjected  to  the  required  amount  of  compressive  load / 
unload  strain  to  induce  the  required  residual  stress  profile  at  each 
particular  depth.  The  simulated  inclusion  sizes  ( <  10  pm)  were 
small  compared  to  the  overall  3D  FEM  model  dimensions,  so  the 
gradient  in  applied  stress  and  residual  stress  over  the  inclusion 
was  considered  to  be  negligible  in  their  analysis.  Alternatively, 
inclusion  sizes  in  powder  metallurgy  polycrystalline  Ni  base 
superalloy  IN100  are  on  the  order  of  10  100  pm  [6,23,24],  For 
inclusions  of  these  sizes,  the  gradient  in  residual  stress  (RS)  field 
due  to  shot  peening  (Ref.  RS  profile  in  [25])  can  have  a  significant 
effect  on  stress/strain  response  and  should  be  considered  when 
analyzing  inclusion  and  RS  effects  in  Ni  base  superalloys.  There 
fore,  a  simulation  approach  that  accounts  for  the  entire  distribu 
tion  of  residual  stress,  and  not  just  at  discrete  surface  depths,  is 
warranted.  Hence,  we  aim  to  develop  a  framework  to  assess 
(1 )  the  effect  of  microstructure  on  the  entire  residual  stress  profile 
due  to  shot  peening  and  (2)  the  effect  of  cyclic  loading  on  residual 
stress  relaxation  in  polycrystalline  Ni  base  superalloy  components. 

This  section  begins  with  a  brief  overview  of  previous  methods 
to  impose  residual  stresses  within  components.  Next,  the  eigen 
strain  application  of  residual  stresses  within  a  FEM  framework  is 
discussed  and  a  J2  plasticity  model  is  presented  for  calibrating  this 
model.  Finally,  the  method  by  which  crystal  plasticity  is  incorpo 
rated  is  presented  with  results  for  variability  of  initial  residual 
stress  and  retained  residual  stress  due  to  a  single  load/unload 
sequence. 


3.2.  Previous  methods  for  simulated  application  of  residual  stresses 

Techniques  to  simulate  the  shot  peening  process  can  be  divided 
into  two  methods:  (1)  Simulation  of  the  impact  response  between 
the  shot  bead  and  the  shot  peened  surface  by  quasi  static  or 
explicit  dynamic  analyses  to  predict  the  resulting  residual  stress 


Material  Pre  Shot  Peening 


Material  -  Post  Shot  Peening 


Idealized  Average  Deformation  of 
Material  Post  Shot  Peening 


Schematic  of  shot  peening  process  Idealization  of  shot  peening  process 

(multiple  impingements).  (single  deformation  event). 


Fig.  1.  Schematic  showing  idealization  of  shot  peening  process  as  a  single 
deformation  event. 


and/or  Almen  intensity  as  a  function  of  shot  peening  parameters 
(shot  size,  speed,  coverage,  material  properties,  incident  angles, 
etc.)  [26  42]  and  (2)  Simulation  of  the  overall  induced  mechanical 
response  due  to  shot  peening  through  a  deformation  process 
[20,43,44],  Some  examples  of  relevant  works  regarding  these 
two  methods  are  discussed  below. 

3.2.1.  Simulation  of  single  and  multiple  impact  events 

The  first  means  to  model  the  high  velocity  impact  shot  peening 
process  focused  on  single  impact  events  on  an  elastic  plastic  target 
substrate.  Chen  and  Hutchinson  [27,28]  and  Boyce  et  al.  [26]  found 
that  explicit  dynamic  simulations  incorporating  effects  of  strain 
rate  sensitivity,  inertia,  and  elastic  wave  propagation  resulted  in  a 
better  prediction  of  residual  stresses  than  quasi  static  analyses. 
Frija  et  al.  [29]  used  an  energy  equivalence  expression  and  3D 
quasi  static  FEM  analysis  to  find  good  correlation  between  the 
predicted  residual  stress  along  the  shot  bead  impact  centerline 
and  experiments.  Zion  and  Johnson  [34]  studied  the  impact  of  a 
hard  and  soft  shot  bead  and  a  high  strength  steel  and  found  that 
the  most  highly  significant  input  factor  was  the  value  of  friction 
coefficient  assumed  between  the  shot  bead  and  target  material. 
While  these  single  impact  simulations  can  unveil  key  relationships 
between  shot  peen  input  (e.g.  shot  diameter,  impact  velocity, 
incident  angle)  and  residual  stress  output  [32],  multiple  impact 
events  should  be  used  to  simulate  more  realistic  peening 
conditions. 

One  means  to  model  multiple  impacts  is  to  combine  discrete 
element  modeling  (DEM)  and  FEM  to  determine  RS  profile.  DEM 
simulates  spatial  interactions/collisions  among  multiple  discrete 
particles  within  the  shot  stream.  DEM  exports  particle/substrate 
impact  velocity,  location,  and  contact  forces  into  a  FEM  model  to 
calculate  resulting  plastic  strains  and  residual  stresses.  Multiple 
researchers  [39  42]  have  used  this  combined  DEM  FEM  approach 
to  link  the  effect  of  complex  part  geometry  on  the  overall  shot 
peening  process,  including  the  effects  of  impact  angle,  impact 
density,  and  coverage  on  location  specific  residual  stress 
formation. 

These  combined  DEM  FEM  simulations  of  multiple  impinge 
ment  events  are  certainly  noteworthy.  However,  it  is  hard  to 
discern  how  much  variability  in  local  residual  stress  is  attributed 
to  (1)  the  localized  inhomogeneous  deformation  due  to  the 
stochastic  impingement  events  inherent  in  the  shot  peening 
process  and  (2)  the  effect  of  local  microstructure,  especially  grain 
orientation.  Therefore,  to  isolate  the  effect  of  local  microstructure, 
we  propose  to  use  a  uniform  average  “amount  of  impingement” 
within  a  crystal  plasticity  finite  element  model.  In  this  approach, 
the  individual  shot  peen  events  that  were  modeled  in  the  previous 
method(s)  are  not  modeled.  Instead,  the  resulting  material  defer 
mation  and  hardening  states  due  to  shot  peening  are  (1)  induced 
explicitly  within  an  implicit  or  explicit  FEM  model  and  (2)  used  as 
initial  conditions  for  subsequent  relaxation/fatigue  analysis 
simulations. 

3.2.2.  Techniques  to  induce  overall  mechanical  response  due  to  shot 
peening 

One  way  to  induce  residual  stresses  explicitly  within  an  FEM 
model  is  to  deform  the  FEM  model  in  displacement  controlled 
constrained  compression.  This  method  mimics  the  collective 
mechanical  means  in  which  biaxial  residual  stresses  are  imposed 
within  a  component  during  the  shot  peening  process  and  traces 
the  evolution  of  the  material  state  throughout  the  deformation 
process.  For  example,  Prasannavenkatesan  et  al.  [20,21]  used  a 
simple  displacement  controlled  method  in  conjunction  with  iso 
tropic  plasticity  [20]  and  polycrystal  plasticity  [21  to  induce 
residual  stresses  [20]  and  reproduce  experimental  trends  in 
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residual  stress  relaxation  [21]  due  to  HCF  cyclic  bending  of  a 
martensitic  gear  steel.  For  these  analyses,  individual  3D  FEM 
models  were  used  at  discrete  depths  and  subjected  to  the  required 
amount  of  compressive  load/unload  strain  to  induce  the  required 
residual  stress  profile  at  each  particular  depth.  This  displacement 
controlled  method  presented  by  Refs.  [20,21]  is  best  used  for 
simulated  application  of  residual  stresses  to  the  surface  of  a 
smooth  planar  specimen  with  well  defined  boundary  and  loading 
conditions. 

An  alternative  means  to  initialize  residual  stresses  within  a 
FEM  model  is  to  use  the  residual  stress  and  material  hardening 
states  for  target  initial  conditions.  For  example,  Buchanan  et  al. 
[25,45]  initialized  their  material  model  with  an  initial  material 
hardening  state  (effective  plastic  strain,  backstress)  and  residual 
stress  state  based  on  experimental  XRD  and  cold  work  measure 
ments  at  different  depths  in  a  Ni  base  superalloy  IN100  material 
[45],  Similarly,  Benedetti  et  al.  [46,47]  imposed  initial  local  yield 
strength  and  hardening  variables  at  different  depths  based  on 
the  experimental  microhardness  measurement  at  that  depth. 
Additionally,  the  initial  residual  stress  profile  was  induced  by 
introducing  an  experimentally  fit  eigenstrain  distribution 
within  the  FEM  model.  The  method  of  introducing  eigenstrains 
within  an  FEM  model  to  produce  residual  stresses  is  covered  in 
more  detail  next. 

J .2.3.  Eigenstrain  method  of  imposing  residual  stresses 

Eigenstrains  have  long  been  used  to  model  micromechanical 
misfit  strains  that  are  not  associated  with  globally/externally 
applied  mechanical  loading.  Some  examples  where  eigenstrains 
are  used  include  internal  stresses  due  to  inclusions/particles  or 
fibers  [48],  differences  in  coefficient  of  thermal  expansion  of 
different  phases/layups  [49],  phase  transformations  [50],  and  heat 
treatment  effects  [51],  Several  authors  have  used  the  eigenstrain 
method  to  model/reconstruct  residual  stresses  induced  by  shot 
peening  [44,52  54],  laser  shock  peening  [55  61],  and  welding 
[62  66].  For  shot  peening  analysis,  the  amount  of  residual  stress 
induced  as  a  function  of  depth  can  be  controlled  by  specifying 
spatial  distributions  of  thermal  eigenstrains.  The  actual  residual 
stress  application  process  is  nominally  isothermal,  so  the  applied 
temperature  change  and  thermal  expansion  eigenstrains  are 
fictitious;  they  are  merely  introduced  as  a  means  to  induce 
residual  stresses  within  a  component  under  isothermal  conditions. 

The  main  challenge  of  this  approach  is  determining  the  correct 
eigenstrain  distribution  required  to  produce  a  given  residual  stress 
profile.  General  solutions/frameworks  to  the  so  called  “inverse 
eigenstrain  problem”  have  been  developed/presented  by  many 
authors  [44,52,53,67  69],  Universally,  these  approaches  assume 
that  the  eigenstrain  distribution  can  be  reconstructed  as  a  super 
position  of  a  truncated  series  of  basis  functions,  i.e„ 

N 

£*(X)  =  ]T  Cj  Zi(X)  (1) 

i  1 

In  this  equation,  N  is  the  number  of  basis  functions,  ^(x),  and  c,- 
are  the  coefficient  multipliers  for  each  basis  function.  The  benefit 
of  this  model  is  that  one  is  at  liberty  to  choose  the  total  number  of 
basis  functions  and  the  form  of  each  basis  function;  as  a  result, 
there  are  multiple  sets  of  basis  functions  that  can  describe  a  given 
eigenstrain  distribution.  For  example,  one  could  use  a  series  of 
smooth  basis  functions  [67  69],  multiple  polynomials  [44],  or 
even  a  superposition  of  multiple  kernel  density  functions  includ 
ing  triangular  functions  [53]  or  normal  distribution  functions  [52] 
as  the  overall  eigenstrain  distribution  estimator.  Regardless  of  the 
functional  form,  the  most  important  factor  that  should  be  con 
sidered  is  how  to  solve  for  the  basis  function  coefficients  and 


whether  the  given  eigenstrain  distribution  is  able  to  reconstruct 
the  desired  residual  stress  profile.  For  example,  Korsunsky  [44] 
used  axisymmetric  plate  theory  to  analytically  find  stresses  and 
deformations  arising  due  to  peening  and  found  the  necessary 
eigenstrain  as  a  function  of  plate  depth  using  polynomial  functions 
as  basis  functions. 

Methods  that  incorporate  both  the  dynamic  material  response 
during  residual  stress  application  and  the  subsequent  effect  that 
these  residual  stresses  have  on  material  behavior  typically  require 
separate  analyses  due  to  different  time  scales  (dynamic  shot/laser 
peen  process  versus  quasi  static  LCF/HCF/creep  loading).  For 
example,  Achintha  et  al.  [55  57]  combined  explicit  dynamic  and 
implicit  FEM  analyses  to  model  the  laser  shock  peening  in  a  Ti 
6A1  4  V  alloy  using  an  eigenstrain  approach  and  an  assumed 
elastic  perfectly  plastic  material  model.  They  found  that  although 
eigenstrain  distributions  were  similar  for  different  specimen 
thicknesses,  the  resulting  residual  stresses  were  quite  different 
among  different  specimen  thicknesses. 


1.3.  Objectives,  scope  and  limitations  of  current  work 

The  previous  eigenstrain  methods  mentioned  in  the  foregoing 
all  suffer  from  the  same  shortcomings:  they  used  simple  elastic 
plastic  models  that  are  unable  to  determine  the  effect  of  local/ 
random  microstructure  (e.g.  grain  size/orientation)  on  residual 
stress  profile  variability  and  residual  stress  relaxation.  Flence,  the 
objective  of  this  paper  is  to  develop  a  combined  eigenstrain  and 
crystal  plasticity  finite  element  simulation  approach  to  address 
this  critical  need. 

There  are  many  approximations  made  with  regard  to  the  actual 
shot  peening  process  in  our  simulations.  The  process  of  shot 
peening  involves  multiple  random  shot  indentations  inducing 
surface  roughness  (cf.  [15,47,70]).  As  stated  previously,  we  do 
not  model  individual  shot  peening  events.  Rather,  the  collective 
shot  peening  process  is  modeled  as  a  uniform,  quasi  static 
displacement  event  (cf.  Fig.  1)  that  spatially  varies  as  a  function 
of  surface  depth.  The  variability  in  residual  stress  comes  into  play 
due  to  the  microstructure,  not  due  to  incomplete  coverage  or 
random/sporadic  shot  locations.  In  practice,  the  XRD  residual 
stress  measurements  are  averaged  over  the  irradiated  area  deter 
mined  by  the  size  of  the  X  ray  beam.  The  experimental  XRD 
average  residual  stress  data  used  for  comparison  in  this  study  was 
found  over  an  irradiated  area  of  3  mm  x  5  mm  [45],  which 
encompasses  over  10,000  grains.  Therefore,  the  highly  localized 
variations  in  residual  stresses  due  to  surface  undulations  are 
averaged  out  for  the  reported  XRD  residual  stress  data.  Therefore, 
we  do  not  simulate  the  random  surface  undulations  due  to  the 
multiple  surface  impingements  induced  during  the  shot  peening 
process.  Instead,  we  impose  a  uniform  “amount  of  impingement” 
in  each  surface  layer  by  means  of  an  applied  eigenstrain. 

Free  surface  effects  with  respect  to  in  plane  stress  components 
were  neglected/avoided  by  assuming  that  the  finite  element 
material  is  from  the  center  of  the  specimen.  The  increased  hard 
ening  of  the  surface  layer  of  the  shot  peened  material  is  accounted 
for  in  the  crystal  plasticity  finite  element  model  via  an  increase  in 
dislocation  density  (and  the  subsequent  increase  in  yield  strength 
[71])  at  the  surface. 

It  should  also  be  noted  that  the  intense  amount  of  near  surface 
cold  working  and  dense  dislocation  network  produced  during  shot 
peening  can  invoke  near  surface  plasticity  induced  refinement  of 
the  microstructure  [45,70],  This  microstructure  refinement  is  not 
accounted  for  in  this  work.  Thus,  the  purpose  of  our  model  is  to 
simulate  the  overall  induced  mechanical  response  due  to  shot 
peening,  rather  than  individual  impact  events  or  grain  refinement. 
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2.  Methodology  for  imposing  residual  stresses  using 
eigenstrain  approach 

The  methodology  used  to  impose  residual  stresses  within  an 
eigenstrain  framework  and  its  extension  to  the  crystal  plasticity 
finite  element  method  are  covered  in  this  section.  The  quasi 
thermal  expansion  distribution  can  be  optimized  to  fit  any 
targeted/measured  residual  stress  profile.  The  actual  process  is 
isothermal,  so  the  thermal  expansion  fields  are  simply  used  to 
induce  plastic  strain  gradients  and  residual  stresses.  Additionally, 
this  technique  has  the  added  benefit  that  it  can  be  used  for  more 
complex  material  response  (e.g.,  crystal  plasticity)  and  more 
complex  geometries  (e.g.,  notches). 

2.1.  Material  model 

2.1.1.  Experiments  used  to  calibrate  model 

During  the  shot  peening  process,  equibiaxial  compressive 
residual  stresses  are  introduced  near  the  surface  due  to  con 
strained  plastic  deformation.  To  satisfy  equilibrium,  tensile  stres 
ses  form  within  the  subsurface  of  the  material.  These  tensile 
stresses  may  decay  with  increasing  depth  for  thicker  specimens 
[72  74]  or  approach  a  steady  state  value  within  the  bulk  of  the 
material  for  thinner  specimens  [25,45],  For  example,  Buchanan 
et  al.  [25]  performed  residual  stress  relaxation  studies  of  2  mm 
thick  powder  metallurgy  (PM)  Ni  base  superalloy  IN100  (25  pm 
average  grain  size)  specimens  that  were  shot  peened  to  an  Almen 
intensity  of  6  A.  The  residual  stress  profile  given  in  [25]  contains 
average  residual  stress  values  measured  over  a  3  mm  x  5  mm 
irradiated  X  ray  region  at  given  depths  from  the  surface.  These 
measured  average  residual  stress  values  were  fit  to  a  curve  of  the 
form  [45] 

0*s(X)=[(ffs  <?int)  +  CiX]exp(  C2X)  +  c 7in[  (2) 

where  the  least  square  values  of  as=  879.0,  (Tjnt  =  205.7, 
Ci  =  67028,  and  C2  =  20.89  were  determined  by  Buchanan  [45] 

to  fit  the  experimental  residual  stress  data  as  a  function  of  depth 
(x)  from  the  surface.  Eq.  (2)  assumes  that  the  residual  stress 
approaches  a  steady  state  internal  value  of  oiM  with  increasing 
depth  (x)  and  is  applicable  for  depth  values  from  the  surface  (x=0) 
to  half  depth  (x=l  mm),  where  the  condition  of  half  symmetry  is 
assumed. 

Since  the  work  of  Buchanan  [45]  contains  information  on  both 
the  initial  residual  stress  curve  and  relaxation  of  residual  stresses 
with  fatigue  loading  and  is  applied  to  thinner  (2  mm)  specimens, 
we  will  use  Eq.  (2)  for  the  initial  residual  stress  curve.  Additionally, 
we  limit  our  computational  specimen  thickness  to  2  mm  to  ensure 
that  it  is  consistent  with  experiments. 

2.1.2.  Isotropic  J2  plasticity  model 

A  rate  independent  J2  plasticity  finite  element  model  with 
combined  isotropic/kinematic  hardening  was  used  to  calibrate 
the  quasi  thermal  residual  stress  application  approach.  The  result 
ing  calibrated  eigenstrain  distribution  was  later  used  as  an  input 
for  a  crystal  plasticity  finite  element  method  (CPFEM)  model.  The 
J2  plasticity  model  employed  is  an  existing  ABAQUS  [75]  material 
model  that  employs  the  Von  Mises  yield  surface 
F  =f(a  a]  rrys  =  0,  where  F=  0  during  plastic  flow,  a  and  a  are 
the  stress  and  back  stress  tensors,  respectively,  and  nys  is  the  Von 
Mises  equivalent  yield  strength.  The  function  /  is  given  by 

f(a  a)  =  yj(S  a)  :  (S  a)  where  S  is  the  deviatoric  stress 
tensor.  An  associated  flow  rule  is  assumed  as  £-  =  gpJE  Here,  tp  is 

the  plastic  strain  rate  tensor  and  ^  =  w'|ep  :  ep  is  the  equivalent 
plastic  strain  rate.  Cyclic  hardening  of  the  yield  stress  is  accounted 


for  by  the  evolution  of  isotropic  hardening  via 
oys  =  a0+Ks[  1  exp(  tep)]  [76,77],  where  a0  is  the  yield  stress  at 
zero  plastic  strain,  ks  is  the  maximum  change  in  the  size  of  the 
yield  surface,  and  b  defines  the  rate  at  which  cyclic  hardening 
occurs.  The  evolution  of  the  back  stress  tensor  is  characterized  by 
“)  rk«,/  [20,75], 

For  this  study,  we  employ  l<= 2  backstress  terms  (cj=  280,900, 
c2  =  10,178,  r]  =  1163.8,  r2  =  55.65)  to  facilitate  better  stress  strain 
fitting  at  lower  (k=  1 )  and  higher  (k= 2)  strains.  This  computational 
model  was  matched  to  experimental  stress  strain  data  [78]  for  a 
coarse  grained  IN100  Ni  base  superalloy  microstructure.  The  result 
ing  stress  strain  plot  comparing  the  computational  and  experimen 
tal  data  is  illustrated  in  Fig.  2.  The  optimized  isotropic/kinematic 
hardening  parameters  for  this  J2  plasticity  model  are  listed  in  the 
upper  left  hand  corner  inset  of  the  Figure.  As  shown  is  this  Figure, 
this  J2  plasticity  model  mimics  the  cyclic  stress  strain  behavior  well 
for  a  very  complex  loading  history.  Thus,  this  J2  plasticity  model  is 
deemed  sufficient  for  the  calibration  of  the  thermal  expansion 
eigenstrain  method  to  impose  residual  stresses  in  smooth  specimen 
components.  These  results  are  presented  later. 

2.1.3.  Polycrystal  plasticity  framework  with  quasi  thermal  expansion 
eigenstrain 

This  section  describes  how  the  concept  of  quasi  thermal 
expansion  eigenstrain  finite  element  method  (covered  in  the  next 
section)  was  extended  in  the  context  of  polycrystal  plasticity  for 
purposes  of  imposing  a  target  subsurface  residual  stress  field  due 
to  shot  peening.  The  benefit  of  crystal  plasticity  relative  to  the  J2 
plasticity  model  is  that  it  can  characterize  microstructure  varia 
bility;  a  number  of  different  statistically  representative  micro 
structure  instantiations  can  be  simulated  to  address  the 
probabilistic  fatigue  strain  life  distribution.  The  goal  of  this  study 
was  to  develop  a  framework  to  induce  a  full  residual  stress  profile 
within  a  crystal  plasticity  finite  element  framework  to  account  for 
this  microstructure  variability.  Such  a  framework  can  be  used  to 
assess  the  effectiveness  of  shot  peening  in  suppressing  near 
surface  crack  initiation  from  inclusions  located  near  to  the  surface 
of  smooth  specimens  [79]. 

To  incorporate  thermal  expansion  eigenstrain  within  the  crystal 
plasticity  finite  element  method  under  ostensibly  isothermal 
loading  conditions,  a  “quasi”  thermal  expansion  contribution  must 
be  included  within  the  deformation  gradient.  As  shown  in  Fig.  3, 
the  total  deformation  gradient  is  assumed  to  be  multiplicatively 


Fig.  2.  Comparison  of  J2  plasticity  model  in  ABAQUS  to  experimental  data  of  a 
coarse  grain  IN100  cycled  at  650  C.  Experimental  data  are  from  Ref.  [78], 
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decomposed  via 

F  =  Fe  ■  Fp  ■  Ffl  (3) 

where  Fe  denotes  elastic  distortion  and  rigid  body  rotation  of  the 
crystal  lattice,  Fp  accounts  for  plastic  deformation  through  dis 
location  glide  along  crystallographic  planes,  and  Fp  constitutes 
thermal  expansion.  It  is  noted  that  temperature  change  is  intro 
duced  only  to  produce  an  eigenstrain  field  to  induce  initial  residual 
stresses  from  associated  elastic  plastic  deformation.  No  other 
properties  need  be  temperature  dependent.  Isotropic,  linearized 
thermal  expansion  is  assumed,  i.e. 

F®=  V/1+2«A(9  I  (4) 

where  a  is  the  isotropic  thermal  expansion  coefficient,  Ad  is  the 
change  in  temperature,  and  I  is  the  second  rank  identity  tensor. 
The  form  of  Eqn.  (4)  was  chosen  so  that  the  linearized  form  of  the 
Green  finite  strain  tensor  due  to  thermal  expansion,  Es,  is  [80] 


2.2.  Finite  element  model  imposition  of  eigenstrain  field 

To  simulate  the  material  stress  state  after  the  shot  peening 
process,  a  finite  element  model  is  used  and  calibrated  to  the 
experimental  residual  stress  XRD  data  reported  by  Buchanan  et  al. 
[14,25,45,83,84]  on  a  supersolvus  PM  Ni  base  superalloy  IN100 
with  an  average  grain  size  of  25  pm  that  was  shot  peened  to  a 
Almen  intensity  of  6  A.  A  quasi  thermal  method  of  application  of 
residual  stresses  is  chosen  in  this  work  because  it  is  easy  to  apply 
and  calibrate  and  it  is  easily  extended  to  more  complex  models  for 
material  response  (e.g.,  crystal  plasticity)  and  more  complex 
geometries  (e.g.,  notches). 


2.2.1.  Application  to  J2  plasticity  model 

The  general  methodology  for  quasi  thermal  application  of 
residual  stresses  to  a  smooth  specimen  is  shown  in  Figs.  4  and  5 
and  can  be  summarized  as  follows: 


Ep  =  1([Fp]7'fp  l)=aA0I  (5) 

The  rest  of  the  microstructure  sensitive  crystal  plasticity  equations 
for  modeling  deformation  behavior  of  a  coarse  grain  IN100  at 
650  C  follow  those  outlined  by  Przybyla  and  McDowell  [70],  A 
detailed  description  of  the  CPFEM  model  is  omitted  in  this  paper 
to  maintain  brevity.  We  only  highlight  the  incorporation  of  a 
quasi  thermal  expansion  deformation  gradient  (F®)  in  this  section 
to  distinguish  what  is  new/unique  in  the  current  CPFEM  model 
relative  to  the  previous  IN100  model.  More  details  on  the  devel 
opment  of  the  1N100  model  and  the  numerical  implementation 
technique  in  ABAQUS  can  be  found  in  the  work  of  Shenoy  et  al. 
[78,81]  and  McGinty  [82],  respectively. 


Intermediate  Thermal  Intermediate  Stress 


Configuration  Free  Configuration 

Fig.  3.  Multiplicative  decomposition  of  the  deformation  gradient,  including  quasi- 
thermal  expansion. 


1.  Initial  configuration:  Fig.  4  shows  the  finite  element  model  used 
to  simulate  residual  stress  application  to  a  smooth  specimen. 
The  experimental  smooth  specimen  [45]  had  a  nominal  gage 
section  length,  width,  and  thickness  ofy=20mm,  z=10mm 
and  x=2  mm,  respectively  (see  Fig.  4  for  x,  y,  and  z  directions). 
Similar  to  the  finite  element  model  employed  by  Buchanan 
et  al.  [25,45],  a  small  portion  in  the  center  of  this  gage  section 
was  used  for  the  finite  element  model.  Half  symmetry  within 
the  depth  (x  dimension)  of  the  material  was  employed  so  that 
the  overall  dimensions  of  the  finite  element  model  for  J2 
plasticity  simulations  were  xdim  =  lmm,  ydim= 34  pm,  and 
zdim= 34  pm.  For  eigenstrain  distribution  calibration,  the  finite 
element  model  was  divided  into  many  2.5  pm  thick  finite 
elements,  as  shown  in  the  rightmost  image  in  Fig.  4.  It  should 
be  noted  that  since  the  compressive  region  of  residual  stress 
field  is  very  thin  (x  depth  ~  200  pm)  and  there  is  a  high 
gradient  of  residual  stress  with  depth,  a  very  fine  finite  element 
thickness  of  5  pm  is  required  near  the  surface  to  provide 
convergence  of  the  FEM  response.  Buchanan  [45]  reported  a 
similar  FEM  mesh  size  requirement  for  convergence.  In  this 
work,  a  finer  mesh  of  2.5  pm  was  used  to  provide  more 
eigenstrain  distribution  data  points  (black  dots  in  Fig.  5)  to 
improve  functional  form  fitting;  this  fitting  is  covered  later. 
Each  finite  element  was  assigned  a  given  quasi  thermal  expan 
sion  coefficient  aj  and  J2  plasticity  material  properties.  The 
initial  distribution  of  aj  values  was  assigned  so  that  the 
resulting  residual  stress  values  were  within  +  250  MPa  from 
the  target  residual  stress  profile  at  a  given  depth  to  avoid  any 
numerical  instabilities  in  the  numerical  optimization  scheme 
described  below. 

2.  Eigenstrain  application:  with  all  surfaces  constrained  from 
normal  displacement,  an  eigenstrain  (.e*hermj  =  ajAT )  is  intro 
duced  within  the  model  as  shown  in  the  upper  left  hand  corner 


Fig.  4.  Finite  element  model  used  to  apply  residual  stress  to  a  smooth  specimen. 
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Fig.  5.  Methodology  for  quasi-thermal  application  of  residual  stresses  to  a  smooth  specimen. 


of  Fig.  5.  It  should  be  noted  that  an  arbitrary  value  of  AT  =  1 
was  applied  uniformly  throughout  the  whole  specimen,  with 
out  loss  of  generality,  since  temperature  is  not  a  physical  field 
variable  during  the  imposition  of  residual  stress. 

3.  Release  surface  constraint:  the  surface  constraint  at  x=0  is 
removed  to  simulate  the  spring  back  of  the  material  after  shot 
bead  impact  (cf.  lower  left  hand  portion  of  Fig.  4).  This  step  is 
required  to  relax  the  stress  component  normal  to  the  surface  so 
that  cjxx  ~  0,  which  is  representative  of  the  near  surface  stress 
state  at  the  end  of  the  shot  peening  process  [20].  This  relaxa 
tion  step  is  analogous  to  the  “unload”  step  in  the  displacement 
controlled  method  of  residual  stress  application  presented  by 
Prasannavenkatesan  et  al.  [20], 

4.  Optimization  of  thermal  expansion  coefficients:  a  secant  root 
finding  method  [85]  is  used  in  conjunction  with  the  FEM  model 
to  optimize  the  spatial  distribution  of  thermal  expansion 
parameters  to  fit  the  experimental  residual  stress  profile. 

5.  Fit  thermal  expansion  coefficients  to  functional  form:  Gaussian 
probability  density  functions  and  polynomials  are  used  to  fit 
the  optimized  thermal  expansion  coefficient  as  a  function  of 
depth  from  the  surface.  More  details  on  this  functional  form 
and  how  it  is  used  within  the  crystal  plasticity  framework  are 
covered  in  the  following. 


the  secant  root  finding  method  and  the  red  and  green  solid 
lines  show  the  functional  form  fitting  of  these  computationally 
optimized  thermal  expansion  coefficients.  The  top  middle  plot 
in  Fig.  5  displays  the  entire  distribution  of  thermal  expansion 
coefficient  as  a  function  of  x  distance  from  the  surface.  The 
bottom  right  two  plots  in  Fig.  5  illustrate  zoomed  in  versions  of 
the  fitting  of  the  required  thermal  expansion  coefficient  to  the 
piecewise  functional  forms  described  below.  A  piecewise 
smooth  functional  form  was  required  so  that  the  thermal 
expansion  coefficients  can  be  defined  independent  of  mesh 
size  and  so  that  this  functional  form  could  be  used  as  an  input 
for  the  crystal  plasticity  model.  The  thermal  expansion  func 
tion,  a(x),  was  split  into  4  sections: 


1.  x  <  0.058  mm:  this  section  of  the  curve  was  fit  using  a  super 
position  of  two  Gaussian  probability  density  functions  (PDFs).  A 
similar  description  of  using  two  Gaussian  PDFs  was  demon 
strated  in  [52]  to  describe  the  eigenstrains  induced  by  shot 
peening  a  GW103  magnesium  alloy.  The  functional  form  used 
to  fit  the  thermal  expansion  coefficients  at  x  <  0.058  mm  is 


a(x)  =  Aexp 


(x  KBf 

2(j^ 


+Dex  p 


(X  PE? 

2rjp 


+  Gx+H 


(6) 


2.2.2.  Fitting  thermal  expansion  coefficient  to  functional  form 

The  right  hand  side  of  Fig.  5  shows  the  required  thermal 
expansion  as  a  function  of  x  distance  (depth)  from  the  surface 
that  is  used  to  replicate  the  target  residual  stress  profile  shown 
in  the  upper  right  hand  portion  of  Fig.  5.  The  black  dots  in  Fig.  5 
show  the  optimized  thermal  expansion  coefficient  found  using 


2.  0.058  mm  <  x  <  0.082  mm:  this  section  was  fit  using  a  5th 
order  polynomial: 


a(x)  =  a21x5  +a22x4  +  a23x3  +  a24x2  +  a25x+ a26  (7) 
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3.  0.082  mm  <  x  <  0.5  mm:  this  section  was  fit  using  a  10th  order 
polynomial  with  coefficients  ranging  from  highest  to  lowest 
ranked  polynomial  terms  as  a31,a32,  •••  ,a39,a310. 

4.  x  >  0.5  mm:  The  last  section  was  described  by  a  constant 
thermal  expansion  coefficient,  a41 . 

it  should  be  stated  that  the  selection  of  the  x  value  bounds  on 
these  different  curves  were  selected  with  the  constraint  that  each 
section  pieced  together  would  represent  the  overall  optimized 
thermal  expansion  distribution.  Each  piecewise  function  mini 
mized  the  Euclidean  norm  (sum  of  squares  error)  within  its  x 
value  bounds.  The  constants  of  the  two  Gaussian  PDFs  curve  were 
found  using  a  Gauss  Newton  numerical  approach  [85]  and  the 
constants  of  the  polynomial  functions  were  found  using  the  built 
in  MATLAB  “polyfit”  function  [86],  The  constants  used  for  each 
curve  section  are  listed  in  Table  1  and  the  resulting  fitting,  again,  is 
shown  on  the  right  hand  side  of  Fig.  5.  These  curve  demarcations 
and  functional  form  bounds  were  chosen  for  convenience.  Any 
number  of  functional  forms  could  be  used.  For  example,  one  could 
use  Eureqa  [87,88]  to  give  any  number  of  solutions  to  this 
eigenstrain  distribution.  The  main  purpose  of  doing  these  demar 
cations  is  so  that  we  can  input  these  functional  forms  into  the 


Table  1 

Constants  used  to  fit  functional  form  for  thermal  expansion  coefficients  as  a 
function  of  specimen  depth. 


Curve  1 

Curve  2 

Curve  3 

Curve  4 

Var. 

Value 

Var. 

Value 

Var. 

Value 

Var. 

Value 

A 

0.00051 

a21 

-4021918 

031 

-5450 

O41 

-0.00227 

HB 

0.0241 

a22 

147844 

032 

14866 

Oc 

0.00764 

023 

-26104 

033 

- 17603 

D 

0.00227 

a2  4 

15814 

034 

11861 

He 

0.0377 

«25 

-580 

035 

-5007 

of 

0.0130 

026 

8.58 

036 

1373 

G 

0.0557 

037 

-24.8 

H 

0.0025 

fl38 

27.4 

O39 

-1.79 

0310 

0.0549 

ciystal  plasticity  User  MATerial  (UMAT)  subroutine  and  induce  a 
given  amount  of  eigenstrain  as  a  function  of  depth  from  the 
surface. 

2.2.3.  Application  to  polycrystal  plasticity  model 

Using  a  similar  finite  element  model  as  previously  described  for 
the  J2  plasticity  model,  a  combined  ciystal  plasticity  and  J2  plasticity 
model  was  constructed  as  illustrated  in  Fig.  6.  The  polycrystalline 
grain  structure  within  the  crystal  plasticity  region  is  constructed 
using  a  random  sequential  adsorption  algorithm  similar  to  that 
described  in  [71,89,90].  This  spherical  packing  algorithm  offers  more 
control  over  grain  sizing  as  compared  to  a  traditional  random  seed 
Voronoi  tessellation,  which  results  in  a  normal  distribution.  The 
values  of  p=  0.1  and  n=0.4  were  chosen  for  the  target  lognormal 
grain  size  (mean  grain  size =34  pm)  distribution  function, 

/(x;  p, a)  =  x^2iexp ^  (lni]gV'>2j .  based  on  previous  publications  of 
fine  grain  IN100  grain  size  distributions  [71,90  92],  An  example  of 
the  target  grain  size  distribution  and  the  actual  grain  size  distribution 
created  using  the  spherical  packing  algorithm  is  shown  in  Fig.  6(a). 
These  grain  size  distributions  are  normalized  by  the  mean  grain 
volume,  ( Vgm )  =  4/3^(0.034  mm)3  =  1.65  x  10~3  mm3. 

Fig.  6(b)  shows  an  isometric  view  of  an  example  polycrystalline 
grain  structure  used  for  the  FEM  application  of  residual  stresses.  In 
this  Figure,  each  grain  is  represented  by  a  different  color  to 
visualize  the  grain  structure.  The  FEM  model  is  a  square  prism  of 
material  with  a  cross  section  that  is  0.17  mm  by  0.17  mm,  which 
corresponds  to  having  approximately  five  grains  through  the  y 
and  z  thicknesses  of  the  cross  section.  Crystal  plasticity  is 
employed  for  elements  that  are  within  0.35  mm  of  the  surface 
and  J2  plasticity  is  employed  for  elements  that  are  at  a  distance 
greater  than  0.35  mm  from  the  surface  to  the  total  x  dimension, 
which  is  1  mm  in  this  case. 

Fig.  6(c)  shows  a  scaled  side  view  of  the  example  polyciystal 
line  grain  structure  shown  in  Fig.  6(b)  overlaid  on  top  of  the  target 
residual  stress  profile  to  compare  the  assigned  FEM  material 
behavior  and  mesh  refinement  to  the  target  residual  stress  profile 
as  a  function  of  depth  from  the  surface.  As  shown  in  Fig.  6(c),  the 
refinement  in  mesh  was  selected  to  correspond  to  key  areas  in  the 


Fig.  6.  Example  combined  polycrystals  plasticity  and  J2  plasticity  finite  element  model  used  for  eigenstrain-based  application  of  residual  stresses.  Crystal  plasticity  is  used  for 
depths  of  x  <  0.35  mm  and  J2  plasticity  for  x  >  0.35  mm.  The  experimental  XRD  residual  stress  profile  is  from  Buchanan  et  al.  [25]  and  the  target  residual  stress  profile  is 
given  by  Eq.  (2).  (a)  Target  versus  actual  gain  size  distribution  using  random  sequential  adsorption  algorithm,  (b)  Isometric  view  of  example  polycrystalline  grain  srtucture 
used  for  FEM  application  of  residual  stress,  and  (c)  scaled  side  view  of  example  polycrystalline  grain  structure  in  (b)  overlaid  on  target  RS  profile  to  show  assigned  material 
behaviour  and  and  mesh  refinement  as  a  function  of  depth  from  the  surface. 
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residual  stress  profile.  Since  a  mesh  size  of  5  pm  provided 
convergence  for  the  FEM  response,  this  mesh  size  was  used  for 
elements  that  were  within  the  compressive  residual  stress  zone, 
x  <  0.2  mm  from  the  surface.  Since  the  experimental  residual 
stress  was  given  up  to  a  depth  of  x=0.35  mm,  the  crystal  plasticity 
zone  was  extended  up  to  this  depth  also.  Beyond  this  depth  of 
x= 0.35  mm,  the  residual  stress  field  is  relatively  flat;  thus,  J2 
plasticity  was  used  for  the  zone  of  x  >  0.35  mm  and  the  mesh  size 
slowly  coarsened  out  to  the  depth  of  x=l  mm. 

2.3.  Method  to  impose  single  load/unload  sequence 

The  boundary  and  loading  conditions  used  during  the  single, 
uniaxial  load/unload  sequence  to  study  residual  stress  relaxation 
are  shown  in  Fig.  7.  Fig.  7(b)  shows  2D  projections  of  the  3D  FEM 
model  (Fig.  7(a))  as  viewed  from  the  surface  (left)  and  side  (right). 
The  loading  step  starts  from  the  condition  where  the  x=0  surface 
constraint  has  been  released  and  the  stresses  in  the  model  have 
been  allowed  to  relax  normal  to  the  surface  (ref.  bottom  left  hand 
side  of  Fig.  5).  Periodic  boundary  conditions  are  applied  to  the 
lateral  surfaces  (at  z=0  and  z=0.17  mm)  using  a  multi  point 
* Equation  constraint  in  ABAQUS  [75],  The  bottom  surface  aty=0 
and  the  surface  at  full  FEM  depth  (x=  1  mm)  are  constrained  from 
normal  displacement  during  the  load/unload  sequence.  Rigid  body 
rotation  and  translation  are  prevented  by  constraining  z  direction 
displacement  for  4  nodes  located  at  coordinates  (x,y,z)  =  {(0,  0,  0), 
(0,  0.17  mm,  0),  (1  mm,  0,  0),  (1  mm,  0.17  mm,  0)};  these  4  node 
locations  are  denoted  in  Fig.  7(b)  by  the  green  filled  black  circles. 
With  these  boundary  conditions  applied,  the  top  surface 
(y=0.17  mm)  of  the  model  is  subjected  to  a  given  net  normal 
traction,  (ryynet=oa ,  and  then  unloaded  to  oyy  net  =  0  MPa.  During 
this  loading  process  the  top  surface  is  subjected  to  a  multi  point 
constraint  (MPC)  to  make  all  of  the  nodes  on  the  top  surface 
displace  the  same  throughout  the  deformation  process. 

A  MPC  displacement  controlled  method  is  used  for  residual 
stress  relaxation  studies  because  residual  stress  relaxation  is  most 
pertinent  to  strain  controlled  conditions,  e.g.,  notches.  The  mag 
nitude  of  residual  stress  relaxation  will  vary  depending  on  the 


maximum  load  or  displacement  applied.  Consequently,  we  study 
the  effect  of  maximum  load/displacement  on  residual  stress 
relaxation. 


3.  Results  and  discussion 

3.3.  Contour  plots  of  stress  in  specimens 

Fig.  8(b)  (e)  shows  example  stress  contour  plot  results  for  the 
polycrystalline  grain  structure  depicted  in  Fig.  8(a).  The  left  hand 
column  of  this  Figure  shows  contours  of  the  rryy  component  of 
stress  and  Von  Mises  stress  at  the  end  of  eigenstrain  application 
with  all  finite  element  surfaces  constrained  (upper  left  of  Fig.  5) 
and  the  right  hand  column  shows  these  same  contour  plots  after 
the  surface  constraint  is  released  so  the  outer  tractions  become 
zero  (lower  left  of  Fig.  5).  For  comparison  purposes,  the  contour 
plot  scales  of  the  oyy  component  plots  (Fig.  8(b)  and  (c))  are 
identical  as  well  as  the  contour  plot  scales  for  the  Von  Mises  plots 
(Fig.  8(d)  and  (e)).  These  Figures  display  the  ability  of  the  FEM 
model  to  induce  compressive  residual  stresses  using  constrained 
thermal  expansion. 

The  first  thing  that  is  immediately  noticed  from  this  contour 
plot  is  the  high  value  of  ayy  stress  (or  more  precisely,  the 
hydrostatic  stress)  induced  near  the  surface  as  the  eigenstrain  is 
applied  and  the  faces  are  constrained  (Fig.  8(b)).  This  high  value  of 
hydrostatic  stress  is  typical  for  a  component  that  is  loaded  in 
constrained  compression.  It  is  also  well  known  that  when  a 
component  is  loaded  in  constrained  compression,  larger  strain 
(or  equivalently  stress)  is  required  for  yielding  as  opposed  to  an 
unconstrained  compression  condition.  In  fact,  constrained  com 
pression  can  increase  the  apparent  yield  strength  by  a  factor  of 
two  as  compared  to  unconstrained  compression.  In  these  con 
strained  compression  cases,  a  better  indicator  of  plastic  response  is 
to  use  a  deviatoric  (or  equivalent)  stress  measurement.  Hence, 
both  the  ayy  and  Von  Mises  stress  measures  are  shown  in  these 
plots.  Although  during  the  constrained  compression  step  the 
maximum  compressive  stress  is  ayy=  3868  MPa,  the  Von  Mises 


Fig.  7.  Schematic  of  boundary  and  loading  conditions  used  for  single  load/unload  residual  stress  relaxation  studies,  (a)  3D  FEM  model  showing  polycrystalline  grain 
structure,  and  (b)  2D  projection  of  FEM  model  as  viewed  from  the  surface  (left)  and  side  (right). 
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are  in  units  of  MPa. 


Fig.  8.  Example  stress  contour  plot  results  for  conditions  before  and  after  surface  constraint  released.  All  stresses  are  in  units  of  MPa.  (a)  Initial  polycrystalline  grain  structure 
(0.17  mm  x  0.17  mm  square  cross-section),  (b)  Contour  of  ayy  at  end  of  eigenstrain  application  with  all  faces  constrained,  (c)  Contour  of  <TyyafLer  surface  constraint  has  been 
released,  (d)  Contour  of  VM  stress  at  end  of  eigenstrain  application  with  all  faces  constrained,  and  (e)  Contour  of  VM  stress  after  surface  constraint  has  been  released. 


equivalent  stress  is  at  a  more  reasonable  maximum  value  of 
1875  MPa. 

The  release  of  the  surface  constraint  on  the  x=0  face  (ref.  the 
bottom  left  portion  of  Fig.  5)  allows  the  finite  element  model  to 
expand  in  the  negative  x  direction.  Subsequently,  the  stress 
component  in  the  x  direction  (cr**)  tends  toward  zero  and  the  ayy 
and  azz  components  settle  into  the  desired  biaxial  residual  stress 
values.  It  should  also  be  noted  that  after  relaxation  of  the  normal 
surface  traction  on  x=0  face,  all  slices  normal  to  the  surface  also 
have  net  zero  ct**  traction  into  the  depth. 

3.2.  Scatter  in  residual  stress  among  multiple  realizations 

In  this  section,  the  results  from  N=  5  instantiations  are  pre 
sented  to  determine  the  amount  of  scatter  exhibited  among 


multiple  digitally  created  microstructures.  Due  to  the  different 
orientation  distribution  of  the  grains,  the  residual  stress  value  of 
an  element  at  a  given  depth  will  differ  depending  on  the  orienta 
tion  of  the  grain  in  which  the  element  is  located  and  its  interac 
tions  between  its  neighbor  grains.  Therefore,  certain  combinations 
of  microstructures  will  cause  the  overall  residual  stress  profile  to 
deviate  above  or  below  a  given  ideal/target  mean  stress  value. 
Plotted  in  Fig.  9(a)  are  results  from  one  random  microstructure 
instantiations.  In  this  Figure,  each  red  dot  indicates  the  residual 
stress  value  within  a  single  element.  Clearly,  there  is  significant 
scatter  in  the  elemental  residual  stress  values,  especially  in  the 
near  surface  area  where  compressive  residual  stresses  are  highest. 
Since  the  finite  element  model  employs  a  structured,  voxelated 
mesh,  the  residual  stress  component  values  at  each  finite  depth  (x 
distance)  are  averaged  and  depicted  with  a  blue  and  green 
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a  b 


Fig.  9.  Computational  versus  experimental  residual  stress  profile  as  a  function  of  depth  from  the  surface.  Experimental  data  are  from  Buchanan  et  al.  [25  ].  (For  interpretation 
of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.)  (a)  Instantiation  1.  and  (b)  Scatter  in  RS  profile  for  5  computational 
instantiations  versus  scatter  in  experimental  RS  data. 
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Fig.  10.  Simulated  residual  stress  relaxation  due  to  a  single  load/unload  sequence  as  a  function  of  maximum  applied  stress.  The  experimental  residual  stress  target  relaxation 
data  measured  after  one  cycle  are  from  Buchanan  [451(a)  Axial  direction  (oyy).  (b)  Transverse  direction  (oy). 


triangle.  Additionally,  the  target  residual  stress  profile  from 
Buchanan  [45)  is  depicted  by  the  thick  black  error  bar  lines.  The 
error  bar  lines  indicate  the  scatter  in  average  residual  stress  from 
six  different  specimens,  where  the  average  residual  stress  from 
each  specimen  is  found  over  an  irradiated  area  of  3  mm  x  5  mm 
[45).  It  is  clear  from  this  plot  that  the  eigenstrain  method  of 
computationally  inducing  residual  stresses  within  the  crystal 
plasticity  framework  results  in  residual  stress  profiles  that  are 
able  to  match  experimental  mean  values. 

Fig.  9(b)  compares  the  experimental  data  from  six  baseline  XRD 
measurements  [45 J  to  the  statistical  spread  from  N= 5  random 
instantiations  using  the  polycrystal  plasticity  model.  The  typical 
definition  of  the  standard  deviation  from  a  small  population  are 
invoked  for  these  ±  3<t  values.  As  shown  in  Fig.  9(b),  there  is  a 
strong  correlation  between  the  statistical  spread  of  the  mean  value 
found  from  experiments  and  the  statistical  spread  of  the  mean 
value  predicted  by  the  crystal  plasticity  simulations  for  the  near 
surface  (x<0.1  mm)  residual  stress  values.  Based  on  only  N=  5 
instantiations,  the  computational  model  slightly  underpredicts  the 
scatter  for  locations  further  into  the  depth  of  the  smooth  speci 
men.  However,  the  current  and  future  studies  focus  on  the  role  of 


residual  stresses  in  suppressing  crack  formation/propagation  near 
the  surface.  Therefore,  we  conclude  that  the  method  used  here 
adequately  supports  these  current  and  future  residual  stress 
studies.  The  utility  of  this  model  will  further  be  explored  with 
respect  to  residual  stress  relaxation  in  the  following  section. 

3.3.  Residual  stress  relaxation  with  single  load/unload  results 

Here  we  compare  the  simulated  relaxation  of  residual  stresses 
to  experimental  residual  stress  relaxation  studies  on  a  coarse 
grain  1N100  at  650  °C  with  average  grain  size  of  25  pm  [45]. 
Buchanan  [45]  presented  the  residual  stress  relaxation  of  two 
specimens  in  the  axial  (y  loading  axis)  direction  and  transverse  (z) 
direction  due  to  a  single  stress  controlled  load/unload  sequence 
with  <rmax  =  900MPa  (Ref.  (45,  Fig.  511).  The  spread  in  XRD 
measurements  of  experimental  residual  stress  relaxation  in  the 
axial  and  transverse  directions  are  used  for  comparison  to  the 
simulated  residual  stress  relaxation  results  presented  here. 

Pictured  in  Fig.  10  is  the  retained  ayy  and  a &  residual  stress 
values  as  a  function  of  depth  and  maximum  applied  stress 
following  the  simulated  load/unload  sequence.  The  top  portion 
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of  this  Figure  shows  the  relaxation  of  the  axial  [ayy)  stress  and  the 
bottom  portion  shows  relaxation  of  the  transverse  (<rzz)  stress  due 
to  the  single  load/unload  sequence.  As  in  previous  Figures,  the 
initial  residual  stress  is  denoted  by  the  thick  black  error  bars.  Also, 
the  gray  shaded  regions  in  Fig.  10  (a)  and  (b)  correspond,  respec 
tively,  to  the  axial  and  transverse  target  residual  stress  relaxation 
profiles  measured  by  XRD  after  unloading  [45], 

The  overall  trend  of  the  relaxation  curves  is  as  expected:  there 
is  an  increase  in  residual  stress  relaxation  as  the  maximum  applied 
stress  is  increased.  When  the  maximum  stress  reaches  a  value  of 
approximately  <Tmax  =  1150  MPa,  the  effect  of  residual  stresses  are 
totally  negated  in  the  axial  direction.  As  the  peak  stress  increases 
to  a  value  of  omax=  1200  MPa,  the  axial  residual  stresses  become 
tensile  at  the  surface  and  go  into  compression  further  into  the 
depth  (x  >  0.4  mm,  which  is  not  depicted  in  these  plots)  of  the 
model.  Although  this  depth  (x=0.4  mm)  at  which  the  computa 
tional  model  transitions  from  tensile  to  compressive  stress  is 
slightly  different  than  the  experimental  one  (x=0.2  mm  [14],  the 
fact  that  the  computational  model  predicts  this  reversal  of  residual 
stresses  is  promising. 

Comparing  the  gray  shaded  target  residual  stress  relaxation 
zone  (from  [45,  Fig.  51])  to  the  computational  relaxation  curves, 
there  is  slight  difference  in  curve  shape  at  depths  in  the  range  of 
x= 0.05  mm  to  x=0.2mm.  The  computational  model  predicts 
more  relaxation  (relatively  speaking)  in  this  region  as  compared 
to  the  experimentally  measured  relaxation  in  this  region.  The 
reason  for  this  difference  is  currently  unknown,  but  potential 
reasons  for  this  discrepancy  are  discussed  in  the  “Limitations” 
section  below. 

The  predicted  residual  stress  relaxation  at  the  surface  of  the 
specimen  in  the  range  of  x  <  0.05  mm  seems  to  follow  the  experi 
mental  trend,  despite  the  error  in  the  range  of  x=0.05mm  to 
x=0.2  mm.  In  the  surface  region,  the  peak  stress  required  to  get  in 
the  range  of  the  experimental  residual  stress  relaxation  is  around 
amax  =  \000  1050  MPa  for  the  axial  direction  and  crmax  =  1100  MPa 
in  the  transverse  direction.  Since  there  were  a  range  of  values  of 
peak  stress  that  resulted  in  residual  relaxation  comparable  to 
experiments,  a  single  peak  stress  value  of  nmax  =  1050  MPa  was 
used  for  assessing  RS  relaxation  behavior  scatter. 

The  same  N=5  instantiations  that  were  used  for  the  residual 
stress  profile  scatter  in  Fig.  9(b)  were  used  to  assess  scatter  of 
residual  stress  relaxation.  In  separate  FEM  simulations,  these  finite 
element  models  were  loaded  up  to  a  maximum  stress  value  of 


c^max  =  1050  MPa  and  then  unloaded  to  zero  stress  along  the  y 
direction.  The  average  residual  stress  profiles  were  obtained  in  the 
axial  and  transverse  directions.  From  these  values,  the  statistical 
spread  of  the  axial  and  transverse  retained  residual  stresses  are 
plotted  in  Fig.  11.  In  this  Figure,  the  left  and  right  columns  contain 
the  retained  axial  and  transverse  residual  stresses,  respectively,  for 
maximum  applied  stress  values  of  1050  MPa.  These  simulated 
values  are  compared  to  the  statistical  spread  from  2  separate 
residual  stress  relaxation  experimental  samples  from  the  experi 
mental  XRD  data  reported  in  Buchanan  [45],  Compared  to  the 
experimental  residual  stress  profiles  [45],  the  scatter  in  residual 
stress  of  the  computational  profiles  seems  to  be  within  a  factor  of 
2  of  the  experimental  profiles.  However,  this  conclusion  is  quite 
preliminary  as  there  were  only  a  limited  number  (N=  2)  of  data 
points  from  the  experiments;  it  is  expected  that  with  more  data 
points,  the  scatter  in  residual  stress  relaxation  data  would 
decrease.  Further  advancement  of  the  residual  stress  relaxation 
model  would  be  possible  with  more  experimental  and  computa 
tional  replicas. 

In  the  current  relaxation  study,  we  only  considered  relaxation 
during  the  first  loading  cycle.  It  is  well  known  that  relaxation  with 
fatigue  cycling  occurs  in  two  stages  [21]:  the  majority  of  residual 
stress  relaxation  occurs  during  the  first  cycle  followed  by  gradual 
relaxation  with  continued  fatigue  cycling.  This  two  stage  relaxa 
tion  process  has  been  reported  for  multiple  materials  in  several 
experimental  [93  95]  and  computational  [21,96,97]  studies  on 
residual  stress  relaxation.  In  the  work  of  Prasannavenkatesan  et  al. 
[20,21],  they  considered  the  effects  of  shot  peening  induced  sur 
face  residual  stresses,  pores,  and  hard  and  soft  primary  inclusions 
in  martensitic  gear  steel  on  nonlocal  fatigue  indicator  parameters 
(FIPs)  and  fatigue  crack  formation  near  the  inclusions.  They 
concluded  that  residual  stress  relaxation  could  only  be  modeled 
using  polycrystal  plasticity  [21], 

3.4.  Limitations  of  the  CPFEM  relaxation  model 

There  are  several  limitations  in  this  model  that  can  cause  errors 
in  the  relaxation  simulations.  First,  the  computational  model 
contains  a  domain  decomposition  of  the  material  into  ciystal 
plasticity  and  J2  plasticity  models.  Although  these  two  models 
were  calibrated  to  the  same  set  of  experimental  data,  the  differ 
ence  in  anisotropic  crystal  plasticity  and  isotropic  J2  plasticity  can 
cause  differences  in  material  stress/strain  behavior  over  the  range 


a  b 


Fig.  11.  Simulated  versus  experimental  scatter  in  residual  stress  relaxation  due  to  a  single  load/unload  sequence  for  maximum  applied  stress  amax  1050  MPa  for  N  5 
instantiations.  The  experimental  residual  stress  scatter  data  are  from  Buchanan  [45].(a)  Axial  (ayy)  relaxation,  amax  1050  MPa,  and  (b)  Transverse  (azz)  relaxation,  amax 
1050  MPa. 
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of  several  grains.  For  example,  Fig.  9(b)  shows  that  the  mean 
residual  stress  behavior  averaged  over  a  bin  size  of  x=5pm, 
y= 0.17  pm,  and  z=0.17  pm  can  vary  up  to  ~200  250  MPa.  This 
scatter  further  increases  when  averaging  over  individual  elements 
(Fig.  9(a))  or  the  scale  of  individual  grains.  On  the  other  hand,  a  J2 
plasticity  FEM  model  smears  the  effect  of  material  anisotropy  and 
predicts  a  single  residual  stress  profile  as  a  function  of  depth  with 
no  scatter.  Additionally,  it  should  be  noted  that  these  material 
models  were  fit  to  a  coarse  grain  IN100  with  a  slightly  larger  grain 
size  (34  pm)  than  the  IN100  microstructure  used  for  the  residual 
stress  experiments  (25  pm). 

Another  potential  source  of  error  could  be  due  to  the  fact  that 
this  model  is  a  quasi  static  representation  of  a  dynamic  shot 
peening  process.  Shot  peening  involves  high  speed  collision  of 
many  shot  beads  against  a  surface.  The  current  model  does  not 
take  into  account  such  dynamic  effects  as  elastic  wave  propaga 
tion,  high  strain  rate  effects,  or  inertia  effects;  Incorporation  of 
these  effects  has  been  shown  to  provide  better  residual  stress 
prediction  compared  to  that  of  a  quasi  static  analysis  [26  28], 
Additionally,  the  intense  amount  of  cold  working  and  resulting 
dense  dislocation  network  produced  at  the  surface  can  invoke 
plasticity  induced  refinement  of  the  microstructure  in  the  surface 
layer  [45,98],  Finer  grain  structure  can  affect  the  yield  point/ 
strength  of  the  material  in  the  thin,  highly  plastically  deformed 
surface  layer;  this  change  in  yield  point  could  significantly  alter 
the  elastic  plastic  response  of  the  whole  model.  This  refinement 
in  the  microstructure  was  not  accounted  for  in  this  work.  How 
ever,  the  purpose  of  the  current  model  was  to  simulate  the  overall 
induced  mechanical  response  due  to  shot  peening,  rather  than 
individual  impact  events  or  grain  refinement. 

Regardless  of  these  limitations,  the  current  framework  is  able 
to  reproduce  essential  aspects  of  the  initial  residual  stress  profile, 
scatter,  and  relaxation.  While  the  total  profile  of  the  residual  stress 
relaxation  curves  do  not  agree  precisely  with  experiments,  the 
simulated  relaxation  trends  near  the  surface  correlate  well  with 
experiments.  Therefore,  the  current  quasi  static  eigenstrain  appli 
cation  of  residual  stresses  is  deemed  reasonable  for  assessing  the 
effect  of  residual  stresses  and  microstructure  on  fatigue  variability. 

4.  Summary 

In  this  work,  a  framework  is  presented  for  imposing  shot 
peened  residual  stresses  using  computational  ciystal  plasticity. 
The  residual  stresses  are  induced  by  a  distribution  of  imposed 
quasi  thermal  expansion  eigenstrains.  This  distribution  was  first  fit 
to  an  experimentally  measured  residual  stress  curve  using  an 
isotropic  J2  plasticity  material  behavior  and  then  applied  within  a 
ciystal  plasticity  finite  element  constitutive  model.  Good  correlation 
between  computational  and  experimental  values  were  obtained  for 
(1)  the  initial  residual  stress  profile,  (2)  the  scatter  in  the  initial 
residual  stress  profile  among  multiple  instantiations,  and  (3)  near 
surface  residual  stress  relaxation  trends  for  a  single  load/unload 
cycle.  This  method  of  coupling  the  effect  of  microstructure  and 
residual  stresses  can  be  used  to  investigate  the  effect  of  certain 
microstructural  features  (such  as  the  largest  grain,  inclusions,  pores, 
etc.)  on  microstructure  sensitive  fatigue  estimation. 
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